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Abstract 



We examine the dynamics of solutions to nonlinear Schrodinger/Gross-Pitaevskii equations 
that arise due to Hamiltonian Hopf (HE) bifurcations-the collision of pairs of eigenvalues on the 
' imaginary axis. To this end, we use inverse scattering to construct localized potentials for this 

model which lead to HH bifurcations in a predictable manner. We perform a formal reduction 
from the partial differential equations (PDE) to a small system of ordinary differential equations 
^ ' (ODE). We show numerically that the behavior of the PDE is well-approximated by that of the 

ODE and that both display Hamiltonian chaos. We analyze the ODE to derive conditions for 
the HH bifurcation and use averaging to explain certain features of the dynamics that we observe 
numerically. 

■ 1 Introduction 

0^ ' In the dynamical systems approach to mathematical physics, an important and physically-motivated 

I approach is to consider the behavior of special solutions: fixed points, periodic orbits, and the like. In 

particular, one often wants to know whether a given solution is stable, i.e. whether it can be destroyed 
by introducing a small perturbation to the initial conditions. Moreover, the study of bifurcations has 
shown that instabilities can, in general, occur in a finite numbers of ways. Oscillatory instabilities in 
Hamiltonian systems arise due to Hamiltonian Hopf (HH) bifurcations and have been seen in a great 
number of analytical and numerical studies, as outlined later in the introduction. The aim of this 
paper is to study in detail the nonlinear dynamics that occur in one such system that arises in various 
''{Jj ' applications as a way to get a handle on this phenomenon in general. 

5^ I The nonlinear Schrodinger/Gross-Pitaevskii equation (NLS/GP) 

idi:^ ^ Hi; ~\ijj\'^'iP: H ^ ^dl + V{x), (1-1) 

is important in mathematical physics in (at least) two main contexts. In nonlinear optics, it arises in 
the paraxial approximation for light propagating in a thin waveguide constructed in a material with 
Kerr nonlinearity [DH]. In a Kerr material, the refractive index of light takes the form n = no+n2(\E\^) 
where E represents the electric field. In particular, the electric field is given by 

E{x,z,t) = 3fi(e^('=^-"*V(2:,C))- 

Here z is the direction of propagation along the waveguide, x the direction transverse, and t is time. 
The waveguide is assumed to be thin in the y direction, and the variation in this direction can safely 
be ignored. The potential V{x) represents the contribution due to the geometry of the waveguide, 
and which we assume to be smooth, negative and exponentially localized. The variable ^ represents 
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a scaled version of the propagation distance z. An effective equation for the envelope t\}{x^ C) can be 
derived by the method of multiple scales, and we assume that the independent and dependent variables 
in this equation can be rescaled to obtain the simple form of (jl.ll) dependent on as few parameters as 
possible. Despite the physical meaning of the variable for the remainder of the paper, we shall call 
this variable t to remind us that it is the independent variable of evolution. 

When the sign on the nonlinear term of (|l.ip is reversed (and (. now genuinely represents time), 
the equation describes the state of a Bose-Einstein condensate, a state of matter achievable at extreme 
low temperatures where atoms lose their individual identities and are described by a common wave- 
function 3 . For equation (jl.ll) to hold, the three-dimensional condensate must be strongly confined 
by a steep potential in the two transverse directions y and z so that it assumes a "cigar" shape. The 
term V{x) then represents a less steeply confining potential in the third spatial dimension. 

In both these systems, a fundamental object of study is the nonlinear bound state, i.e. a localized 
solution to (11.11) of the form 

ip{x,t) = e-*"**(a;). 

A solution consists of ^(a;), a sufHciently rapidly decaying real- valued function, and two real numbers 
n and J\f that satisfy 

= - 




'^^{x)dx = \\^\\l=Af 



The parameter TV > 0, the square of the L^-norm, represents the number of particles of a BEC or the 
total intensity of the light in optics. This solution may be thought of as a nonlinear generalization of 
an eigenfunction of a linear Schrodinger equation, although, of course, the principle of superposition 
fails to apply in this instance. We expect, and find, continuous families of solutions that are indexed 
by the intensity Af. In fact, as A/" — >■ 0, some of these solutions approach, in shape, the eigenpairs of 
the linear system. 

Nonlinear bound states, or standing waves, represent coherent and simple states that might be 
observable in a laboratory experiment. Such bound states may be found numerically, or, for specially 
constructed potentials V{x)^ might be easily computed exactly in the linear limit and approximately 
for A/" 7^ 0. In order for such states to be observable in experiments, they would have to be stable, i.e. 
if a solution to equation (jl.ip is initialized at t = with value close to, but not equal to, a solution to 
system (|1.2p . then it must stay in a neighborhood of that solution for all t > 0. 

Much work, of course, has gone into studying the stability, especially the spectral stability, of 
solutions, i.e. the presence of of unstable modes (corresponding to spectrum with positive real part) 
in the linearization of (|1.1|) about a given solution. In particular, we may think of A/" as a bifurcation 
parameter. Since it is usually the case that there exist families of solutions to (|1.2p continuously 
parameterized by J\f, we may ask for what values of Af the solution is stable. 

The stability of a standing wave is not, however, the whole story. Bifurcation theory dictates that 
there is a relatively small set of scenarios (bifurcation types) that may be observed in the transition 
from stability to instability, and in each of these scenarios certain types of solutions and dynamics may 
be observed. System (|1.1|) is Hamiltonian, and this fact further restricts the types of behaviors that 
can be seen near a bifurcation. 

Several recent studies have focused on the types of bifurcations observable in system (jl.ip and 
related systems and we review a few of them here, in order to motivate the current study. In addition 
to the stability of a solution changing as a parameter is varied, a bifurcation may create new solutions. 
Kirr et al., for example, have demonstrated that solutions to (|1.2|) with a double- well potential 

vj^^\x)^V{x^L) + V{x + L) (1.3) 

undergo a symmetry-breaking bifurcation as the parameter Af is raised from zero [4]. At a critical 
value A/sB, a symmetric solution to equation (II. 2p loses stability and two stable, asymmetric standing 
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wave modes are created. Kapitula, Kevrekids, and Chen |5] have shown that for a triple well potential 
of the form 

(-^) ^ y (-3, _ i) + v{x) + V{x + L). (1.4) 

that these symmetry-breaking bifurcations are replaced by saddle-node bifurcations. In the symmetry 
breaking bifurcations, the new families of standing waves "branch off" of the existing families exactly 
at the location of the bifurcation, while for a saddle-node bifurcation, the new families of solutions are 
not connected to the the existing families at this point. 

Also associated with bifurcations are certain features in the dynamics in a neighborhood of the 
family of solutions. The symmetry-breaking bifurcation studied by Kirr et al. was shown by Marzuola 
and Weinstein to display the dynamics typical of such systems. Below the bifurcation, the ODE 
system has a single-well potential energy,and thus a one-parameter family, of periodic orbits. Above 
the bifurcation, the potential energy has a dual- well shape and thus three topologically distinct families 
of periodic orbits. This manifests itself in a wobbling of the shape of the asymmetric solutions or a 
periodic exchange of energy between the two wells [6] ; see also [7l [8] . 

One particular type of bifurcation that can give rise to much more complicated dynamics is the HH 
bifurcation. While [5] concentrates on enumerating all the standing wave states, they also numerically 
compute the stability of these standing waves, and they do demonstrate a HH bifurcation (figure 6d); 
see also [§]. The HH bifurcation has also been observed in other NLS-related settings. Several studies 
have demonstrated numerically the existence of "Krein collisions" — defined in section [5731 below — in 
discrete wave equations [lOl HIl US and in Bose-Einstein condensates (BEG) fTS l [T6 l fTZl [T8] . 
In these studies, and most others, the bifurcation is discussed only in the context of detecting the 
instability transition in the linear spectrum, or by performing a small number of numerical solutions 
to the initial value problem. 

The HH bifurcation is often described as an instability resulting from a "collision of modes" (i.e. 
frequencies), for example in describing the motion of multiple dark solitons in a quasi-one-dimensional 
BEC's [m fig. 5c], Theocharis et al. remark on an instability caused by "the collision of the second 
anomalous mode with the quadrupole mode" in describing dynamics that look remarkably like our 
figure inSbj column 3. A goal of this paper is to shed light on the origin of such patterns in this and 
similar numerical simulations. 

Kapitula et al. have developed rigorous analytical methods for counting the number eigenvalues 
that might lead to instability in a wide variety of Hamiltonian nonlinear wave equations [19[ I20j , and 
are thus able to rigorously determine the stability of localized solutions of these infinite-dimensional 
Hamiltonian systems In that work, they apply this method to investigate the stability of localized 
solutions to a system of coupled NLS equations. In [3T], they use this machinery to study the stability of 
rotating matter waves in Bose-Einstein condensates, and demonstrate the presence of HH bifurcations. 
They supplement this with well-chosen numerical simulations in order to demonstrate the dynamics 
that occur when the solution is destabilized. 

In related work, Goodman and Weinstein 22 study the linear stability of standing wave modes 
of the nonlinear coupled mode equations (NLCME). In that paper, several bifurcation scenarios are 
outlined, including both symmetry-breaking (figure 4.2c) and HH (figure 4. 2d). In extensive numerical 
studies, they found symmetry-breaking bifurcations, but were unable to locate any HH bifurcations. 
Part of the motivation for the construction in the present paper was to engineer potentials where these 
bifurcations can be observed and understood, first in the simpler and better-known NLS/GP equation. 
In forthcoming work parallel to this, we perform similar analysis for NLCME and find largely similar 
results. 

In this paper, we focus the dynamics in the vicinity of a HH bifurcation. In the following subsection, 
we summarize the notation used in the paper. In section [21 we discuss the assumptions about the 
potential under which this bifurcation may be observed and state the main findings of this paper, 
including a slight reformulation of the problem in section 12.31 In section [3l we sketch the inverse- 
scattering techniques used to construct the potential, while leaving more of the details to Appendix [Xj 
Section [4| discusses the elementary properties of the finite-dimensional model. In section [4TT1 we briefly 
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describe the derivation of a finite-dimensional model 14.61 for the of the dynamics of equation 11.11 and 
in section 14.31 a further reduction 14.141 of the dimension based symmetries of the system. Section 14.21 
reviews the known stationary solutions of svstem lTBl In section[5l we derive a formula to detect the HH 
bifurcation. Section|6]contains numerical confirmation of this formula and numerical explorations of the 
dynamics of both the PDE Il.ll and the finite-dimensional model (|4.14p . We discuss a further symmetry 
reduction of the system (|4.14|) in section [71 which allow a fuller understanding of the dynamics, and 
finally conclude in section [51 Appendix IbI contains some formulas related to the derivation in section [31 

1.1 Notation 

• An overbar, z represents the complex conjugate of z. 

• The expressions 9fiz and represent, respectively, the real and imaginary parts of z. 

• We denote the L? inner product over complex-valued functions of a real argument by (/, g) = 
S^f{x)9{x)dx. 

2 Technical Background 

2.1 Discrete spectrum of the operator H. 

If V{x) has even symmetry, V{—x) — —V{x), then solutions to the (linear) eigenvalue equation 

VL^^H^, (2.1) 

that is, the A/" — ?> limit of equation (|1.2|) . will will have either odd or even symmetry. If the NLS 
system (|2.1|) possesses two discrete eigenvalues Vli < VI2 < 0, then standard Sturm-Liouville theory 
requires that the associated eigenfunctions ^1 and ^'2 are, respectively, even and odd functions of x. 
^'i is the minimizer of the associated Hamiltonian and is thus referred to as the ground state. The 
mode ^'2 is referred to as the excited state. The spectrum of H will, independently of its symmetry, 
generically consist of a finite number of real discrete eigenvalues fife < and continuous spectrum on 
the non-negative real axis. 

These standing wave modes persist as M is increased from zero, with their shapes and frequencies 
altering as well. For sufficiently small amplitudes, they will inherit the neutral stability of their linear 
limits-barring resonances among the the eigenvalues fife that we will discuss shortly. In [4], Kirr et al. 
prove that as the norm of the solution is increased, then at a critical amplitude 

A/'sB oc - r^i, 

the solution that continues from the ground state loses stability and a new stable solution to ()1.2p 
appears, possessing neither even nor odd symmetry. That is, there is a symmetry-breaking or (Hamil- 
tonian) supercritical pitchfork bifurcation. Marzuola and Weinstein have demonstrated for this system 
in the unstable regime, over a long time period, the dynamics of (|l.ip are well- approximated by a Duff- 
ing oscillator-like dynamics [6] when the initial condition is sufficiently close to an elliptic fixed point. 
Pelinovsky and Phan have generalized this result to a wider class of initial conditions and provided a 
proof that relies on simpler estimates [8]. 
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The present problem is naturally modeled by a three degree-of-freedom Hamiltonian system, which 
due to symmetry, as we will discuss, can be reduced to a two degree-of-freedom system. It should also 
be noted that the symmetry-breaking bifurcation is non-generic-if V{x) is non-symmetric, the system 
will generally feature a saddle-node bifurcation instead. HH bifurcations are not possible in the two- 
mode system, and to observe them, we must consider a system with an additional degree of freedom. 
We demonstrate via formal asymptotics, and observe numerically, that 5*2, the first excited state, 
generically becomes unstable in an HH bifurcation when V{x) supports three localized eigenmodes 
and with eigenfrequencies satisfying the following assumptions. 

Assumptions 

(Ai) ni<n2<^3< 0, 

(A2) 172-^^1 = 0(1), 
(A3) 173-1^2 = 0(1), 
(A4) (1^3 - ^2) - {^2 - ni) < 1, and 

(A5) ~ O (1) (i.e. a sufficient gap between the three eigenmodes and the band edge). 

To satisfy assumption | ( A4)] in particular, we let 

1^2 - 111 = - e and - = + e (2.2) 

where e <^W and W = O {1). The sign of e is left unspecified while > 0. 

By using inverse scattering techniques, we can construct a potential V{x) with whatever eigenvalues 
we choose and which also satisfies the evenness condition. In fact, the HH bifurcation is generic and 
will occur regardless of the evenness of V{x). The behavior of the system above the critical amplitude 
may, however, affect the nonlinear behavior of the system in the supercritical regime. 

2.2 Symmetries 

Let m{i;) ^H'tP'- |?Ar V- Then for any real potential V{x), *Tt possesses 0(2) symmetry. More specif- 
ically, defining the operators R^f{x) — e*'^/(a;), and Zf{x) = f{x) corresponding to multiplication by 
an arbitrary complex phase and complex conjugation, we see that 

fn(i?^V) = i?0fn(V') and "RiZtij) = zfn(V'). (2.3) 

Finally, define the operator R^{f{x)) = f{—x). If, in addition, V{x) is an even function, 01 is also 
equivariant to the Z2 operation 

R-m{i^{x)) ^ m{R^ijj{x)). (2.4) 

Putting these together gives shows that system has 0(2) x Z2 symmetry. Bifurcations in systems 
with such symmetries generally have codimension greater than or equal to bifurcations in similar 
systems without such symmetries. Earlier studies have noted that the results can be generalized to a 
larger class of nonlinearities for which OT(^) is equivariant under (|2.3p and (12. 4p . The same is almost 
certainly true in the present case as well. We choose to work with the simple cubic nonlinearity 
described above because a more general nonlinearity would invalidate relation (|4.5p below and increase 
even further the number of terms in equation (|4.6p . As in Kirr et al., we will show that the reduced 
ODE system has the same symmetries. 
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2.3 An alternate formulation 



If we make the change of variables ip — \/M ijj in equation (ll.ip and = \0V ^ in (|1.2p , we get the 
modified evolution equation 

idt^p = HiP -Af\i;\'^4>, (2.5) 

and stationary equations 

/°° (2 6) 

^^{x)dx = 11*11^ ==1. 
-oo 

This formulation presents a natural environment for studying the A/" — > limit. Since this system is 
well-defined regardless of the sign of J\f, we can study all the bifurcations for Af G M., which gives a 
fuller picture of the dynamics, unifying the focusing and defocusing NLS equations. In section |4l we 
derive finite-dimensional models of systems (|l.ip and (ll.2p . A similar change of variables will allow 
us to put a small parameter N of either sign in front of the the nonlinear terms in, for example, 
system (|4.6p and other equations derived from it. Also, it should be noted, that in this formulation 
there will generally be no bifurcation at Af — 0: for almost all potentials V{x), a smooth family of 
functions will pass right through any solution to system (j2.6p with M = 0. 



2.4 The language of stability and resonance 

Suppose that J\f — in the systems (|2.5p and (|2.6p and that the linear eigenvalue problem has n 
linearly independent solutions (\I'„,il„). Ther0 

n 

^(a;,i) = ^c,e-'"^**,(x) (2.7) 

solves equation (j2.5p . In general, this solution is quasiperiodic: each individual component is periodic, 
but in general, the periods will be irreconcilable, and the solution as a whole is non-periodic. Topolog- 
ically, such a solution lies on an rt-dimensional torus T" in the 2n dimensional phase space, which can 
be thought of as the product of n circles or equivalently as an n-dimensional hypercube, with opposite 
(hyper-)faces identified. A resonance relation is a solution to the equation 

n 

kjflj = {k, n) with A? e Z" \ {0}. (2.8) 

The sum 

n 

defines the order of a given resonance. For example under assumption (|2.2p with e = 0, the vector 
k = (1, —2, 1) satisfies equation (|2.8p and defines a resonance of order 4. The number of independent 
solutions of equation (|2.8p with a given order defines the multiplicity of a that resonance at that order. 
If the system has no such resonances, then each solution (12.71) is dense on T". The number of linearly 
independent vectors A that solve equation (j2.8p is the multiplicity of the resonance. If the system is 
resonant with multiplicity to, then the solutions are confined to, and dense on, n — m-dimensional 
subsets of T" which are themselves topologically equivalent to T"^™. To understand this, think of 
the two-dimensional case. If there are two non-resonant frequencies, the solution is dense on a two- 
torus, so its closure is two-dimensional, but if the two frequencies are rationally related, then each 

^In a finite dimensional model, the eigenfunction 'l'fe(x) would be replaced by an eigenvector v^^'K 
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one-dimensional orbit will be closed. These closed orbits must lie on the level set of an additional 
conservation law. When there is a near-resonance, 

n 

but non-zero, and nonlinear terms are nonzero but small, there will be a nearly conserved quantity 
that allows us to use averaging to decrease the dimension of the system and obtain simpler equations 
that are valid for a finite time. 

Solutions of equation (j2.6l) (but not equation (|2.7|) !) are known as relative fixed points. Simply put, 
when viewed in an appropriate reference frame oscillating with frequency fi, they are time- invariant. 
Similarly, there may exist relative periodic orbits^ which are themselves quasi-periodic, but appear 
periodic when viewed in an appropriate reference frame. 

The linear stability of a given solution to some general system will be determined by the eigenvalues 
Xj of a certain matrix M. The imaginary parts of the eigenvalues determine the frequency with 
which small perturbations oscillate about the solution, and the real parts will determine the growth 
(3?Aj > 0) or decay rate (3?Aj < 0) of perturbations. The solution is therefore unstable if there exist 
any eigenvalues Xj with j > 0. Points in parameter space where the stability changes are called 
bifurcation points, and there are different types. The manifestation of the symmetry-breaking (or 
Hamiltonian pitchfork) bifurcation for NLS/GP is discussed in great detail in [3J|S]. We will describe 
the HH bifurcation in greater detail in section [5751 In Hamiltonian systems, it is well-known that if A 
is an eigenvalue, then so are —A, A, and —A. This implies that the eigenvalues can occur in four types 
of groupings, up to multiplicity: complex quadruplets {A, A, —A, —A} with nonzero real and imaginary 
parts, real-valued pairs {A, —A}, purely imaginary pairs {i/i, — i/i} and zero eigenvalues of even algebraic 
multiplicity. The symmetry-breaking, or Hamiltonian pitchfork occurs when, as a parameter is varied, 
a purely imaginary pair of eigenvalues collide at the origin, producing a purely real pair. Here, small 
perturbations to the origin will initially grow monotonically due to the real positive eigenvalue. See 
figure I2.ir a) and (b) . The Hamiltonian Hopf bifurcation occurs when two pairs of pure imaginary 
eigenvalues collide at a nonzero point on the imaginary axis, and the four eigenvalues recombine to 
form a quartet of fully complex eigenvalues. The dynamics in the near-linear regime is oscillatory due 
to the imaginary parts of the eigenvalues; see figure I^TT c) and (d). The major goal of this paper is to 
investigate the behavior of solutions in the nonlinear regime. 

3 Construction of the linear potential V{x) 

By explicit construction following Harrell 24 , Kirr et al. demonstrate that, given a potential V{X) 

(2) ______ 

with exactly n discrete eigenmodes, the dual well potential V£ i^) given by (|1.3|) will have exactly 
2n eigenmodes, and that the eigenvalues come in pairs, each pair exponentially close to each other and 
to the corresponding eigenvalue of the single-well potential. Kapitula et al. discuss this same idea for 
a three- well potential given by ()1.4|) . As L — > oo, the three eigenvalues of this system all converge to 
a multiplicity-three eigenvalue — a highly degenerate situation, for which a complete analysis is rather 
more complicated, with many bifurcations occurring quite near to each other. 

Another way to proceed is to specify the eigenvalues flj = —K^j = l...n and to use inverse 
scattering methods to construct a reflectionless potential with exactly these eigenvalues [25]. This 
will be unique except for n integrating factors £_j (corresponding to the positions of the solitons) that 
arise from solving the associated Gel'fand-Levitan-Marcenko equations. The solution will exactly be 
a two-soliton solution of the Korteweg-de Vries equation [26'. Backlund transformations and Darboux 
transformations can be used to more easily find this solution. The Darboux transformation is very 
similar to the Backlund transformation, but it yields not only the potential, but its eigenvectors, which 
will be useful in what follows [27l[28]. There is a unique way to choose the constants such that V{x) 
is an even function (corresponding to the situation where all n solitons collide at the origin at t = 0). 
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Figure 2.1; (a) The path of the eigenvalues as an parameter is varied in the Hamihonian pitchfork 
bifurcation, (b) The real and imaginary parts of the eigenvalues, (c) The path of the eigenvalues as an 
parameter is varied in the HH bifurcation, (d) The real and imaginary parts of the eigenvalues. After 
Luzzatto-Fegiz and Williamson ^23) . 



When n — 2, the general formula for this 2-soliton is 

4(k2 — k^)(k2 cosh2Kia; + cosh2Koa;) 
((ki — K2) cosh (k2 + Ki)x + {ki + K2) cosh (ki — K2)x) 

with Ki > K2 > 0. This has (un-normalized) ground state and excited states 

cosh K2X 



and 5*2 = 



(ki — K2) cosh (ki + K2)x + (ki + K2) cosh (ki — K2)x 
sinh Kix 

(ki — K2) cosh {k2 + Ki)x + (ki + K2) cosh (ki — K2)x 



and frequencies Qj = When ki — 2 and K2 = 1, this potential reduces to the familiar initial 

condition for the KdV two-soliton 

V{x) — —6 sech^ X 



with frequencies fii = —4 and ^2 = —1- If we choose ki = vT+e and K2 — \/l — e, then for 
< e ^ 1, the potential (13.11) takes the form of dual- well potential, very similar to that studied by 
Kirr et al. However the eigenvalues and eigenfunctions are now known exactly. 

To compute the Hopf bifurcation, we may construct a potential V{x) as a three-soliton solution to 
KdV. The three-soliton has a very similar form to the two-soliton in equation (|3.ip , but with ten terms 
in the numerator and four in the denominator. It is also computed via the Darboux transformation. 
This solution is given in appendix [A] If the parameters are chosen such that the three eigenvalues 
are spaced very closely together (close to a triply-degenerate eigenvalue) , then this potential takes the 
form of three nearly identical potentials spaced equidistantly apart at a large distance, as was studied 
by Kapitula et al. In a similar, but much more complex, vein, Hirsh et al. have used inverse scattering 
in order to design potentials that support modes of a user-prescribed shape I29j . 

An example that displays the HH bifurcation, and which we will use in our subsequent numerical 
studies, is shown in figure [5711 Here the potential is chosen with il = (—11.1,— 10,— 9.1). In fact. 
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for the rest of this paper, this potential will be used, except where otherwise noted. It is the mode 
corresponding to — 10 that undergoes the HH bifurcation. 

4 The finite-dimensional model 
4.1 Derivation of the model 

We decompose the solution to equation (|2.5|) as the following time-dependent linear combination: 

= ci(t)^'i(a;) + C2{t)-^2{x) + ca-^aix) + ri{x; t) (4.1) 

where the eigenvectors are orthonormal and, for each t and j, ri{x]t) is in the orthogonal complement 
to the discrete eigenspace, i.e. 

= Sij and (?7(-,t),*j) = 0, for i,j = 1,2,3. 
We define the projection operators on to the discrete eigenmodes 

n,C-(*„C>*., for J- = 1,2, 3 (4.2) 

and onto the continuous spectrum 

HcontC = C ~ (Hi + Ha + n3)C. (4.3) 

Following the methodology of Marzuola and Weinstein, we substitute the decomposition (|4.ip into 
the PDE (|l.ip and apply to it the four projection operators defined above, giving evolution equations 
for the components of the decomposition. The following system of equations is equivalent to the 
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PDE (j2.5p under the assumption that V{x) is even and supports exactly three modes. 

2 2 2 2 2 

aiiii|ci| ci + aiii3(ciC3 + 2 |ci| C3) + aii22(2ci |c2| + C1C2) 

+ aii33(2ci |C3|^ + C1C3) + ai223(c2C3 + 2 |C2|^ C3) + 01333 |C3|^ C3 



QlCi + Af 

at 



i?l(ci,C2,C3,?7) 

(4.4a) 



at 



ail22(c?C2 + 2 |ci| C2) + 2ai223(ciC2C3 + C1C2C3 + C1C2C3) 

+02222 |C2|^ C2 + a2233(2C2 jcsl^ + C2C3 



*$-f^3C3+AA 

at 



fliiialcil ci + aii33(ciC3 + 2 |ci| C3) + ai223(2ci |c2| + C1C2) 

+ ai333(2ci |C3|^ + C1C3) + a2233(c2C3 + 2 |c2|^ C3) + 03333 |c3|^ C3 



i?2 (ci,C2,C3,?7) 

(4.4b) 



^3(ci,C2,C3,77) 

(4.4c) 



IC^tT? - Hl]+M\l]\ rj = i?cont(ci,C2, 03,77); 

(4.4d) 



where 



where we have used that if {ttj, Tr^, tt;, TTm} is any permutation of {j, fc, m}, then 



(4.5) 



and that ajkim = if j + fc + / + to = 1 mod 2. These parameters will be calculated numerically 
as needed for the simulations presented below. The Rj and -Rcont terms are the projections onto the 
appropriate eigenspaces of remaining nonlinear terms of (j2.5l) and are presented in full in appendix [Bl 
Ignoring the contributions of rj{x]t) to the solution, we derive a finite-dimensional approximation 
to (mH) 



- Q.ICI + N 

at 



aiiii|ci| ci + aiii3(ciC3 + 2 |ci| C3) + aii22(2ci |c2| + Cic\ 



+aii33(2ci |C3| + C1C3) + ai223(c2C3 + 2 |C2| C3) + ai333 |C3| C3 



= 



(4.6a) 



i-^-n2C2+N 

at 



ail22(CiC2 + 2 |ci| C2) + 2ai223(ciC2C3 + C1C2C3 + C1C2C3) 

+02222 |C2|^ C2 + a2233(2c2 |c3|^ + C2C3) 



(4.6b) 



- f^3C3 + N 

at 



aiii3|ci| Ci + aii33(c^C3 + 2 |ci| C3) + ai223(2ci |c2| + CiCj) 



+ai333(2ci |C3| + C1C3) + a2233(c2C3 + 2 |C2| C3) + 03333 jcsl C3 =0; (4.6c) 



We have the following slight change of notation in this equation. System (|4.4[) . being equivalent to 
equation (|2.5|) conserves the norm 



|ci| 



|C2| 



|C3| 



This implies that 



|C1|' + |C2|' + |C3|'<1. 



System (|4.6p possesses a finite-dimensional conserved quantity 



|C1|' + |C2|' 



|C3| =1. 
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The conserved quantities in systems (j2.5|) and (|4.6|) are not equivalent, since the contribution of ri{x, t) 
is ignored in the latter. Recall that Af represents the total intensity and the sign of the nonlinearity 
in equation (|l.ll) . Since the meaning of A/" is slightly changed from equation (|2.5p to system (|4.6|) . we 
introduce the new constant N. 

Note that system (|4.6p possess the same 0(2) x Z2 symmetry as the PDE (jl.ip from which they are 
derived. The 0{2) symmetry is defined in the obvious way, analogous to that used in equation ()2.3p . 
The Z2 symmetry is due to the equivariance of equation (j4.6p generated by the operation of 

^^flip(ci,C2,C3) = (ci,-C2,C3). 

Note that composing this with the operator i?^ with </> = tt flips the signs on ci and C3, leaving C2 
unchanged. 

4.2 Stationary Solutions 

We first look for stationary solutions of system ()4.6p of the form 





This calculation is well-covered by Kapitula et al. [5: in the case where the aijki coefficients satisfy 
some properties that significantly simplify all the equations and the resulting analysis. We will repeat 
the parts we will need in what follows. It is simple to show that the solution to equation ()1.2[) with 
real potential V{x) is, up to a constant phase factor, a real-valued function. This allows us, without 
loss of generality to assume x,y, z € M and that the stationary solution is of the forrcH 

{n - Q,i)x + N{aiiiix^ + 3aiii3X^z + 3aii22xy^ -I- San^sxz"^ + 8012232/^^ + 01333^^) = (4.7a) 

(f2 - fl2)y + N{3aii22x'^ + 6ai223xz + 02222?/^ + 3a2233^^)y = (4.7b) 

{Q, - fi3)z + A^(aiii3a;^ -I- 801133x^2 + 8012232;?/^ + 8ai333a;z^ + 8022332/^^ + 03333^^) = (4.7c) 

x^ + z^-l = (4.7d) 



Odd solutions 

Note that and '03 lie in the invariant subspace of even functions, and thus that this system has 
solutions with y = and x and z nonzero. Similarly, since only -02 fies in the odd invariant subspace, 
there are solutions with only y nonzero. It is these whose stability we will investigate. This family of 
solutions, 

(a^odd, yodd, Zodd, ^odd) = (0, 1, 0, ^2 - 02222^), (4.8) 
may lose stability in a HH bifurcation. 

Even solutions 

li y — then equation (|4.7dp allows us to write x = cos 9 and z — sm9. Eliminating ft from the 
remaining equations yields a single equation 

n(^- ai333 sin'^6' -t- (03333 - 801133) sin^ 6* cos 6* -I- 8(01333 - 01113) sin^6'cos^ 9 

+ (801133 — oiiii) sin^cos^ 9 + 01113 cos^ 9^ + (fli — fl^) sin 6* cos 6* — (4.9) 

■^See Kirr et al. for a more thorough justification [4J. 
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which has period tt, allowing us to consider the domain < 9 < n. For iV ^ 1, this has two solutions: 
one near 6 — and one near = 7r/2 which correspond to the fixed points {x,y,z,U) = (1,0, 0, 
and (x, y, z, 51) = (0, 0, 1, ils) of the linear problem. Both these solutions are found numerically to be 
stable for all N, which is consistent with the Krein signatures of their linearizations. Additional even 
solutions may appear due to saddle-node bifurcations as discussed in [Sj and are shown in figure |4?T] 

General Solutions 

Finally, there are solutions with all three components nonzero. Cancelling out a factor of y in equa- 
tion (j4.7bp . then system (|4.7|) depends on y only through j/^. We therefore define polar coordinates 

X = rcosO, z — rsinO, < r < 1, 

and use equation (j4.7dp to get y^ = 1 ~ r^. Note that the cases r = and r = 1 correspond to even 
and odd solutions discussed above. We use equation (j4.7bl) to solve for in terms of r and 9. After 
plugging this value into the x- and z-equations, we have a system in r and 9 alone. Eliminating r from 
this system produces a single equation for 9: 



where 



= /34 sin'' 9 + 133 sin^ 6i cos 6* -I- /^z sin^ 9 cos^ 9 + pi sin 9 cos^ 9 + CDs'* 9 



(4.10) 



Pa — Af( — 1133312222 + 9ai223a2233 + 3ai333a2233 — 3oi223a3333) — (3ai223 — ai333)(f^2 — ^s)', 

Pz = ^(18a^223 — 9^2233 ~ 6ai223ai333 " 3aii33a2222 + 9aii22a2233 + 9aii33a2233 — 3ii22a3333 + 0222203333)+ 

(02222 — 6a2233 + 03333)^^1 + (— 3aii22 + 3aii33 + 302233 ^ 03333)^^2 + (3aii22 — 3aii33 — 02222 + 302233)^^3 
P2 — A^(27aii22ai223 — 90112201333 — 80111302222 + 80133302222 + 90111302233 — 2701223^2233)+ 

(-9oi223 + 3oi333)fii -|- (801113 - 3oi333)fi2 + ("801113 + 9oi223)f^3; 
Pi = A^(90ii22 — I8O1223 ~ 9oii220ii33 -I- 601113O1223 — 0111102222 + 801133O2222 + 801111O2233 ~ 9oii2202233) + 

(—801122 + 801133I -I- 02222 — 802233)^^1 + (oiiii — 801122 — 3oii33 -1- 802233)^^2 + (— oiiii + 601122 — 02222 

Pq — A^( — 3O1113O1122 + 801111O1223 — 9oii220i223 + O1113O2222) + (01113 — 8oi223)(^^l " ^^2)- 

It is straightforward, though messy to find the saddle-node bifurcation values of N where new 
solutions to equations (|4.9p and (|4.10p arise (although this becomes much neater if the coefficients 
Oijfez are assumed to satisfy condition (8.6) of [Sj. This calculation shows that in both cases, the 
bifurcations happen for = O {W) which we assume to be O (1) as e — t- 0. A complete bifurcation 
diagram of solutions to system (|4.7|) is shown in figure |4?T] This figure does not show more complicated 
solutions to system (|4.6|) such as quasiperiodic orbits whose existence we demonstrate in later sections. 

4.3 Model Reduction via Symmetries 

The HH bifurcation does not lead to any new fixed points. Further, as the Cj may evolve, it is not 
sufficient to assume that each component is real. We therefore must work directly with equation (|4.6p 
rather than with the simpler equation ()4.7p . 

System ()4.6p may be written in Hamiltonian form as 



i/ =lll |ci|V r!2 |C2|' + r!3 |C3|' - A^ 



2'-l'-3 



jOiiiilcil + 01113 |ci| (C1C3 -I- C1C3) 

+ 01122 (\clcl + 2 |ci|' |C2|' + ic?C^) + 01133 (^C^C^ -f 2 |ci|' |C3 

/ 2 2 2 I 2 

+ Oi223 ( 2|c2| (ciC3 + C1C3) + C1C2C3 -I- C1C2C3 1 -f 01333 |C3| (C1C3 -|- C1C3) 

+ ^02222 |C2|** + 02233 (^C2C3 + 2 |c2|^ |C3|^ + 5C2C3) + ^03333 |c3|'' 



(4.11) 



12 




Figure 4.1: The complete set of stationary solutions to system (|4.7p . with the symmetry of the cor- 
responding PDE sohitions indicated by Une style. The amplitudes of the Hopf bifurcations discussed 
in sections [5] and |6] are indicated by points on the curve of odd solutions. The bifurcation creating 
heteroclinic orbits of the even subspace as given by equation (|7.8[) is shown as a horizontal dotted line. 



with evolution equations 

dH 

Because H is equivariant under the group action 

(C1,C2,C3) e'''(ci,C2,C3), 

the dynamics also conserve the (squared) /^-norm: 

|ci|' + |c2p + |c3p = l (4.12) 

as a consequence of the equivariance of system (|4.6|) to the operator given in section 12.21 (see also 
equation (|4.7d|) ). Physically, this simply states that the number of particles or photons is conserved. 
We may use this property to reduce the problem from three degrees of freedom to two. Taking 
advantage of the phase-invariance of H, we define new evolution variables 



Cl 



{t) = ai(Oe*''(*); c^it) = p{t)e''^'^- c^{t) = a3(t)e^^(*). (4.13) 



where p{t),6{t) E M. The Hamiltonian (|4.11l) is independent of 9, and conservation law (I4.12p tells 
us that p{t) = (1 — |cri(/:)| — |o'3(i)| Y^^. Using this we may write down the reduced Hamiltonian 
dependent on just cti, (T3, and their complex conjugates: 

H =(rii - ^2) |criP + {^3 ~ ^2) ksP - N ifliiii IctiI'' + aiii3 |crip (cria-3 + g-iCTs) 

+ aii22(l - Wif - |fT3|')(|fTi|' + 2(3?ai)2) + aiM^a^al + 2 |a3|' + ^afal) 
+ 01223(1 — IfiT — |c3|^)(o'i(T3 + 2aias + 2cti(T3 + d-ia^) + ai333 |o-3|^ {(Jia^ + aias) 

+ 102222(1 - |cti|' - ksl')' + 02233(1 - kil' - |a3|')(|a3|' + 2(3?a3)') + ^03333 ^3! 
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This gives evolution equations 

i&i =(0i - 512)0-1 + N - anil |cri|^cri - 01113(2 \(Ti\'^ as + crfa-a) 

+ 01122 ((2 \aif + |cr3|^ - l)(2(Ti + CTi) + i((T^ - ai)a-i^ - 01133(2 |ct3|^(Ti + ctiCT; 
+ 01223 ((2|ai|^ + 1^31^ - 1)(2(J3 +^3) +CT?(a3 + 2a: 



(4.14a) 



01333 



|cr3|^cr3 - O2222 (Ici |^ + {(T^f - l)cri + ^02233 (4 |cr3 1^ + CT3 +0-3)0-1 



and 



i&3 =(fl3 - fl2)f73 + N - aiii3 IfJil^ (Ti + 501122(4 IctiI^ + crf + af)a3 - 01133(2 [ctiI^CTs + 0-1CT3) 



+ 01223 ((kil' + 2 IfTal' - l)(2ai + cti) + a^((Ti + 2ai)) 
- 01333(2 lasl^ ai + 0-10-3) - O2222(|o-i|^ + |o-3|^ - 1)0-3 
+ 02233 ((ki|' + 2 lasf - l)(2a3 + ^3) + ^(^1 - ^1)^3) 



03333I0-3I 0-3 



(4.14b) 



This reduction involves fixing a reference phase 9(t) and thus leads to equations that are not equivariant 
■with respect to operator R^. 

The full solution may be recovered using the conservation law 



P 



l^W,\'-\a,\\ 



and the auxiliary equation for 0{t): 



m =-n2 + N 02222(1 - Wif ~ Wsf) + ioii22(ki|' + 2^{al)) 



Oi223(2CT3ai + 2ctiCT3 + (Tia3 + CT1CT3) + |02233(|0-3|^ + 2^{a'^)) 



5 Linear stability 



We are particularly interested in the stability of the antisymmetric mode ^2(2;, A^), the nonlinear 
continuation of the linear eigenmode 4*2 of equation (I2.ip . 

5.1 Linearization of PDE solutions 

Letting (^,51) be a solution of system (|2.6p and consider small time-dependent perturbations of the 
form 

ip{x, t) = -^{x) + (m(x, t) + iv(x, t))e^'^K 
Then, linearizing and making the standard assumption that 



u = C/(x)e^S V = V{x)e 



one finds the eigenvalue problem 



A 



V 



n + dl- V{x) + 3JV 



-{n + dl~v{x)+j\f)\ (u 
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5.2 Linearization of ODE 



First, we determine the linear stability of the solution (|4.8[) . which corresponds to ai 
system (|4.14p . By inserting the form 



£73 = in 



ui{t) + ivi{t) 
U3{t) + ivz{t) 



into system (|4.14p . the linearized equations become 



di 



Us 

Vl 
\V3/ 



O2 
M2 



Ml 
O2 



where O2 is a 2 x 2 matrix of zeros 
Ml 

and 



Qi -Q2 + (02222 - 01122)^^ 

-ai223^ 



Vl 
\V3J 



-ai22sN 

^2 + (^2222 — 02233)^ 



(5.1) 



M2 = 

Note in the limit N 



-{VLi - VI2) + (3aii22 - a2222)N 

3ai223^ 



3ai223-/V 

-(fis - i}2) + (3a2233 - a2222)N 



) has eigenvalues ±i{Q2 — ^1) and ±i(il3 — 1^2) on 
the imaginary axis, and that for all N, the matrix M is symplectic, i.e. M = JK where K is symmetric 
and 

'O2 -I\ ^, ^ J. (M2 O2 
/ O2) «°thatX=f ^ 



0, the matrix M 



J 



5.3 Analytical criterion for ODE bifurcation 

The HH bifurcation is the result of a Krein coUsion, where two (pairs of complex-conjugate) eigenvalues 
collide on the imaginary axis, as is shown schematically in figure lOT c) and (d) [30]. If ^ is any vector 
in the eigenspace belonging to the eigenvalues ±jw on the imaginary axis. Then the Krein signature 



m) = sgn 



is constant on the entire eigenspace, and is thus a property of the eigenspace rather than of any 
particular nonzero vector in that space. If two eigenvalues have opposite Krein signature, then, upon 
collision, they will generically split into a Krein quartet, indicating that the origin has become unstable, 
with oscillatory dynamics due to their nonzero imaginary parts. 

In the present ODE, in the — > limit, the eigenvalues, i.e. the frequencies in the linear system, 
are ±i(r22 — ^1) and ±i{il2 — ^3. The Krein signatures are /C(±i(f22 — ^1)) = sgn(i7i — ^2) and 
/C(±j(r23 — r22)) ~ sgn (r^a — 512), implying by assumption ( |(Al)[ ), that their Krein signatures are 
different. Thus, their collision will lead to instability. In fact, the Krein signature can be interpreted 
as the direction of phase rotation, and since 5^2 lies between fli and fis, the Krein signatures can 
be determined without performing this calculation. We expect, and find numerically below, that the 
frequency at which the collision takes place is near ± 1512 ~ ^^3! ~ ± M2 — ^^3! ~ W. 

To detect the HH bifurcation, we construct P{X;N), the characteristic polynomial of M, which, 
as is generic for Hamiltonian systems, is a quadratic polynomial in A^. Letting q = X^, we define 
the simpler quadratic polynomial p{q; N). There will be a double eigenvalue at the value of N where 
the discriminant of p{q) is zero. We further make assumption ()2.2|) . The discriminant is a quartic 
polynomial n(iV). Defining v — N/e, and factoring out a common factor of e^, we find 



n(iy) = di{e)i^^ + dsiey + d2iey + di{e)v + do{t) = 



(5.2) 
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where 

C?4 =(3aii22 — 4a2222 + 3a2233)^(ail22 ^ 20223301122 + 4a^223 + '^2233)^^ 
=8(aii22 — a2233)(ail22 — 02222 + a2233)(3ail22 — 4a2222 + 3a2233)W^e 

— 8(3aii22 — 4a2222 + 3a2233)(ail22 ~ 20223301122 + 4aj223 + 02233)6^ 
d2 =16(aii22 — ai223 ^ 02222 + 02233)(oil22 + 01223 ~ 02222 + 02233)W^^ 

— 8(aii22 — a2233)(7aii22 — 8a2222 + 7a2233)W^e 
+ 16(a^j22 ^ 2a22330ii22 + 4a^223 + 02233)^^ 

di =32(aii22 - a2233)W^e - 32(aii22 - 02222 + 02233)^^^ 

do =\m'^ 

We may solve this numerically or by a perturbation expansion of the form: 

1/ = m + O (e) 
and find that there are double eigenvalues at 

iVHH.± = T + 0r(e2). (5.3) 

—01122 ± 01223 + 02222 — 02233 

Thus, there will be a HH bifurcations for small values of e. The O (e'^) term (not shown) contains a 
factor of W~^, showing that if ^ 1, the divergence of the bifurcation value from a simple linear 
function of e is greater. This is exactly the case of a near-triply degenerate eigenvalue, the case studied 
by Kapitula et al. 

Before proceeding to simulate and analyze system (|4.14p . we make some observations: 

• The first two terms correspond to |cip and |c3p and are unchanged. Similarly any term multi- 
plying a coefficient aijki where each term in the subscript is 1 or 3 is unchanged. Only the terms 
that had a contribution from C2 (those with 2's in the subscripts) are altered. 

• Except for a scaling factor, the real parts of cti and (73 can be interpreted as position variables, 
and their imaginary parts as canonical momentum variables. 

6 Numerical Simulations 

Bifurcation study: spectrum of linearization 

We first consider whether equation (j5.3l) provides a good approximation to the critical value iVun 
when the eigenvalues are of the form given by equation (12. 2p . We show two examples. In both cases 
we choose 1^2 = —10, while e in equation p.2p is allowed to vary. The first subfigure shows W = 1 
and in the second, W = 5. The potential pictured in figure l3T] corresponds to choosing e = 0.1 in the 
first subfigure, and its shape does not change much as e is varied. The result is shown in figure I5TT1 
and demonstrates the large effect W has on higher-order terms in this approximation. We also see 
from this figure, that for small values of e, the system undergoes HH bifurcations at both positive 
and negative values of N, but that for larger values of e, one of these bifurcations may cease to exist. 
This change in character occurs for values of e where two roots of the discriminant ()5.2|) collide and 
annihilate each other. This will happen at values of e where the discriminant of n(7V) vanishes. This 
is formally a sixth-degree polynomial in e, but the coefficients aijki are themselves functions of e, so 
the equation is in fact transcendental. 

Next we investigate how well the bifurcation structure of ODE system (|4.14|) compares with that of 
standing waves of system (|1.2p when the potential V{x) is given in Appendix [A] We find numerically 
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Figure 6.1: Stability diagrams. The trivial solution to system (|5.1|) is linearly stable inside the shaded 
regions, and loses stability to HH bifurcations along the thick black lines. The dashed lines show the 
approximate values (j5.3p of A'hh, computed to O (e). This corresponds to potentials with frequencies 
(a) (r2i,f]2,^^3) = (-11 + 6,-10,-9 + e) and (b) (f^i, 1^2, f^g) = (-15 + e, -10, -5 + e). The vertical 
line in (a) gives the values of the parameters used in all other figures in this paper. 

approximate solutions to (11.21) by first replacing the derivatives with their pseudospectral approxima- 
tions, and solving the resulting (finite dimensional) equations using Matlab's f solve command. We 
then form the linearization about the standing wave, again using the pseudospectral approximation for 
the derivatives. This approximation is implemented as a function and the spectrum is calculated using 
the Matlab eigs command. The discrete spectrum of this problem is compared with the numerically 
calculated spectrum of the matrix in equation (|5.ip . 

When e = 0, the discrete and continuous spectrum of the linearized operator lie entirely on the 
imaginary axis. We can define the Hopf bifurcation point as the value of N for which the eigenvalues 
on the imaginary axis, in pairs, collide and move off into the complex plane as a quartet. In figure [521 
we see that the discrete spectrum of the ODE closely resembles that of the PDE system. The PDE 
system also has a double eigenvalue at the origin, which is unaffected by the bifurcation and which is 
not shown. 




-3 -2 -1 1 2 3 4 5 -4 -2 2 4 6 

N N 



Figure 6.2: The imaginary and real parts of the discrete eigenvalues of the reduced ODE (solid) PDE 
standing wave (dashed). Both show four HH bifurcations. Both the PDE and ODE solutions lose 
stability near = ±0.44 the solutions regain stability N = 4.71 and N = -3.58(PDE) and N = -3.46 
and N = 5.01 (ODE), where it regains stability. The parameter values are ft — (—11.1, —10, —9.1). 
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ODE dynamics 



To simulate system (|4.14l) , we use a Hamiltonian Boundary Value Method (HBVM) of Brugnano et al. 
[311 1321 [33] . which exactly conserves the energy in polynomial Hamiltonian systems. We consider the 
potential shown in figure IXTl which has eigenvalues { — 11.1, —10, —9.1}. The numerically calculated 
value of A^HH = 0.381. We run several numerical experiments, each with initial conditions ai = = 
10~^ and a2 — p — {N — af — cr|), where N S {0.35,0.6,1,3,5.5} corresponding to the five points 
labeled in figure 16.21 These simulations are shown in the five rows of figure 16.31 The first column 
contains a time series of SRcti. The second contains a Poincare section, defined on level sets of the 
amplitude N and the Hamiltonian H, by 

^N,H = {(ai,a3)|9(ai) = & ^5(^1) > 0}. 

Because the values of N, H, and on this section uniquely determine ai S M, T,j^_h is parame- 
terized by the value of a^. We denote the mapping from one crossing oiTiN.H to the next as M. The 
third column contains a reconstruction of the field amplitude computed from ()4.ip and ()4.13p 

In case (a), N — 0.35, the solution oscillates quasiperiodically in a neighborhood of the initial 
condition, with the amplitude of oscillation depending on, and remaining close to, the initial condition. 
A Poincare section, as defined above, shows these solutions appear to be quasiperiodic; see column two 
of the figure. In case (b), N = 0.6, the solution makes large excursions from the initial condition, lying 
close to an apparent homoclinic orbit. In the reconstructed field, column 3, we see that these bursts 
consist of an oscillatory growth of the field in the middle of the potential. We will discuss the exact 
nature of this solution in section[S] For case (c), N — 1.0, the solution appears "weakly chaotic," in 
that it still takes the form of large heteroclinic bursts, only now the time between bursts is irregular. 
The chaos is evident from the Poincare section. In case (d), = 2.25, the chaos is fully developed, and 
the trajectory in Poincare section appears to cover a large open set. This open set avoids, however, 
an elliptical region toward the left half of the figure, on which the map has a fixed point of elliptical 
type. Finally, in case (e) with N — 5.5, we see that for N sufficiently large, the trivial solution is again 
stable. 

At the bifurcation amplitude iVuH, a new fixed point dp of the Poincare map M appears at a 
distance O {\^N — Nhh) from zero, corresponding to a new periodic orbit of the two-degree-of-freedom 
system ()4.14|) . i.e. a new relative periodic orbit of system of the full three-degree-of- freedom system 
defined by the Hamiltonian (|4.11|) . The fixed point Cp appears on the symmetry axis 173 € M. A complete 
periodic solution to reduced system (|4.14l) with N — 2, and a reconstruction of a PDF solution from 
this ODF solution are shown in figure [6^ In subfigure (a), we see that when cti and as are purely real, 
they are in phase. From figure we see that in this case, the modes ^'i and ^3 add constructively 
in the middle well and destructively on the two outer wells. When the ai and (T3 are purely imaginary, 
they are 180° out of phase, so that ^1 and ^3 add destructively in the middle well and constructively 
on the two outer wells. 

As the system reaches the second HH bifurcation a.t N ^ 4.71 where the trivial solution regains 
stability, the new periodic orbit does not disappear, nor does the chaotic motion shown in row (d) of 
figure [6731 Instead, a small region (in fact, a KAM island) around the origin appears at this amplitude, 
on which the solution is regular (quasiperiodic and confined to topological ellipses), and this region 
grows as N is further increased. This is confirmed by numerical simulation. Johansson finds similar 
Hamiltonian chaos when the parameters in his NLS trimer are in the unstable domain, as well as KAM 
islands [1^ 

PDE dynamics 

For comparison, we compute time-dependent solutions of the PDE system. For this we use a Matlab 
code written by T. Dohnal. It uses fourth-order centered differences to compute spatial derivatives, 
and an implicit-explicit additive Runge-Kutta method for time stepping |34j and most importantly for 
long-term simulation, uses perfectly matched layers (PML) to handle the outgoing radiation |35) . 



18 



(b) 



(e) 



100 150 200 250 



10D 200 



400 500 


























0.2 


























0.15 












1 














1 


0.1 
0.05 






'"ill 


'','1 






ill 






|.„i 








l' 








fl 


r")iii'|| 






li' 








llji f" 
























II 




,j -0.05 




























-0.1 




























-0.15 


























;i 



0.05 0.1 0.15 0.2 


0.25 0.3 








\ 




} 


% 


.f 







150 £00 250 300 350 




01. 1 0. Ill Of o;'^i o.s o.;<:i 4 0.45 



200 300 



400 500 



^ll)il|l(l|#IHINW 



-0.4 -0.2 



50 100 




50 100 



Figure 6.3: Simulations of ODE system ()4.14p . The rows, labeled (a)-(e) correspond to the values 
of N indicated on figure |621 The first column shows diai{t). The second shows the intersection of 
the solution with Poincare section Y,n_h- The third column shows a reconstruction of \u(t)\ using 
ansatz (I4.13P (darker areas indicated larger values). This shows, as is predicted by figure \6l2\ that 
cases (a) and (e) are stable, and that instabilities, and even chaos, exist in the other three cases (chaos, 
in fact exists in the fifth case as well). Note in rows (b) and (c) that the Poincare map A4 was run to 
t — 10000 and t — 5000 respectively. Also note the elliptical region toward the left in this figure, into 
which no points of the trajectory enter. For all simulations, the initial condition is ai = CT3 = 10^'^. 
The pictures for nearby initial conditions are qualitatively the same. 
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t 



Figure 6.4: The periodic orbit of the averaged system with parameters as in figure 16.11 with N — 
2 > iVHH- (a) cri(t) and crsit). At <: = 0, the orbits start at the points marked •, and proceed at 
quarter-periods through the points marked ■, A, and ■A'. A reconstruction of the PDE field over three 
periods of oscillation: (b) absolute value, (c) real part (mod e'^*^*^), (d) imaginary part (mod e'^*^*^). 
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As initial conditions, we use linear combinations of the three linear modes, and, to compare the 
solutions with those of the ODE system, we compute the projection of the solution onto the span 
of the localized linear modes, giving us, essentially, the parameters Cj{t). Dividing Cj by the phase 
of C2{t) gives a value analogous to aj and p{t). PDE simulations are shown in figure [6751 again for 
Af G {0.35,0.6,1,3,5.5}. The behavior is remarkably similar to the behavior found for the ODE 
solutions. 

7 Further reduction of the ODE 

We now perform further analysis with the goal of understanding the dynamics observed at amplitudes 
above the critical value for HH bifurcation. We will formally apply the von Zeipel averaging procedure, 
which applies in the case of a resonance between eigenvalues [36i i37j . Because the systems are Hamilto- 
nian, the averaged equations will preserve some but not all features of the full system of equations — for 
example hyperbolic fixed points and their local un/stable manifolds will be preserved, but homoclinic 
orbits will not. The averaged system will be completely integrable, but we have already seen evidence 
in figure l673l that the full system is not. 

The standard reference for the HH bifurcation is the monograph of van der Meer [S^, but the 
analysis presented there does not apply to system ()4.14|) . In the generic HH bifurcation, the matrix 
of the linearization is non-semisimple (i.e. it has non-trivial Jordan blocks) whereas in this case, it is 
semisimple (diagonalizable over C). This particular case is analyzed by Chow and Kim These 
methods are based on a Lyapunov-Schmidt reduction and require a more involved calculation which 
we defer to a later study. 

To reduce the number of degrees of freedom, we first make the change of variables to canonical 
polar coordinates 

c^i VPje'^^-J = 1,3; 

yielding a Hamiltonian: 

-ffpoiar = {N{aii22 COS 26*1 + 2aii22 - 02222) + W - e)pi 

+ {N{a2233 COS 26*3 - 02222 + 202233) - W - e) ps 

+ 2ai223A^(2 cos (6*1 - 6*3) + cos {61 + 6*3)) VP^^/p^ 
+ yp?(-2aii22 cos 26*1 + anil - 4aii22 + 02222) 

- Npip3{aii22 cos 2^1 - aii33 cos 2{9i - ^3) + 02233 cos 26*3 + 2aii22 - 201133 ~ 02222 + 202233) 

+ yP3(-202233 COS26'3 + a2222 - 4a2233 + 03333) 

+ 2A^VP3/0?^^((01113 - 2ai223) COS (6*1 - 6*3) - ai223 COS (6*1 + 6*3)) 

- 2N p^^"^ ^fp{(2ax2i3 COS (6*1 - 6*3) + 01223 cos {Qx + Q3) - 01333 cos (6*1 - 6*3)). 

(7.1) 

Naively, one would hope to make near-identity changes of variables that have the effect of averaging out 
all of the mean-zero (i.e. cosine) terms. Note this would also eliminate the terms of fractional power 
in the pj. The formal equations necessary to remove some of these terms, however, will in some cases 
lead to zero denominators, that is, those terms are resonant. Were it not for such resonances between 
the eigenvalues, one could make near-identity changes of variables to remove all terms containing 
trigonometric functions and fractional powers of p\ and p3 , putting the system in the so-called Birkhoff 
normal form. Resonances of higher order terms become a problem precisely when the linear part of 
the equations contains a resonance of the type defined in equation (j2.8|) . In this case, the linear part 
of the Hamiltonian in the linear limit — )■ 

H = {W - e)pi + {-W - e)p3 = wipi + LJ2P2 
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Figure 6.5: Time-dependent simulations of PDE system (jl.ip . The rows, labeled (a)-(e) correspond to 
the values oi Af indicated on figure Column 1 shows 3?o'i(i). Column 2 shows the intersection of 
the solution with Poincare section Y,j^/^h- Column 3 shows a reconstruction of \u(t)\ using ansatz (|4.13p 
(darker areas indicated larger values). This shows, as is predicted by figure [121 that cases (a) and (e) 
are stable, and that instabilities, and even chaos, exist in the other three cases. The agreement with 
the reduced ODE is uncanny. The initial conditions used are u{x) — O.OOf 'i/'i(a;) + V'2(a^) + 0.001?/'3(a;). 
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satisfies the near resonance of order 2: 

kiUJi + k^Lj^ = 2e ^ 1 

wliere ki = — 1. In tliis case, one cannot completely average the system and is forced to consider 
the Gustavson normal form. For more information see Wiggins [371 §19-10, §20.9]. To put the system 
in normal form, we find (^1,^3) G satisfying kil^ — k^li = 1 and make the symplectic change of 
variables: 

6*1= hi'i-h'^Ps; pi = kiJi+liJ3; 

We also make the assumption that the nonlincarity is small: N — tv. In particular, we could 
choose (^1,^3) — (0,1). This change of variables explicitly separates the fast motion with 0(1) time 
scales from the slower motion with time scales of O (e~^). Such a change of variables would allow us to 
eliminate the pair (Jsjips) from the Hamiltonian. Because in figures l673l and [6?5] we show the Poincare 
map that eliminates the pair (Ji,'0i), we choose to make the equivalent canonical change of variables 

01 = "01, ^3 = -ipi + "03, Pi = Ji + J3, P3 = Js- 
This puts the Hamiltonian in the form 

-Produced = Hq{Ji) + eHi{Ji, J3,^pi,1p3) 

where 

HoiJi) = WJi 

and Hi has period tt in ■01 and 27r in -03 (and which we will not write out here). 

Thus, on a level set of the Hamiltonian iJrcduccd = Wh, we may solve for Ji as a function of the 
other three variables, which gives 

Ji = W + eLi(/i,J3, 01,03) + 0(e2). 

where 



This indicates that to leading order in e and for times of O (e^^), Ji ~ /i is a conserved quantity and 
allows us to use V'l a-s a time-like variable. Renaming ipi = t gives Hamiltonian |40j : 

-^reduced = -eil r^Hi{J3,llj3,h) + ■^7?i( J3, -03, -^l, t) (7.2) 

where 



Hi{J3,^3\ h) = -I1J3+I2JI + 13^/J3\/ J3 + h {2J3 + h~ 1) cos -03 
and Hi, the details of which will not be important, satisfies 



/ ffi(J3,03,-/i,r)dr = 
Jo 



with coefficients 



71 = 2s + {-aiiiih + 2aii22(3ft, - 1) - 201133/1 - 2a2222(^ - 1) + 202233(^1 - 1)) 

72 = (— fliiii + 8aii22 — 4aii33 ~ 4a2222 + 802233 ~ ^3333) 

73 = 2i^ai223- 



23 



Standard averaging techniques [40] now show that there exists a near-identity change of variables 



J= J3 + 0(e),^' = ^3 + 0(e) 
such that the solution to the averaged system with Hamiltonian 



-^average = ^' ^) C^'^) 



agrees with solutions to system (|7.2p with error of order e for times of order . Further, for sufficiently 
small e, fixed points and their local invariant manifolds of system (j7.3p will correspond to periodic orbits 
and their local invariant manifolds of system (|7.2p . 

By the conservation of the total intensity, equation (j4.12p . pi and pa in system (j7.ip are confined 
to the triangle 

< pi < 1; < p2 < 1; < pi + p2 < 1. 
In the reduced system, this becomes a constraint on the conserved parameter h and the variable J, 

1 - h 



-l<h<l; min(-/i,0) < J < 

z, 

We will consider the case Q < h < 1. For the case — 1 < ft, < 0, it is more convenient to eliminate 
(p3, 6*3) and work in the ( Ji, V'l) space. In this case, the phase space is the disk J < 

A short word on this reduction is in order. The level sets of H which are manifolds of dimension 
2n — 1. When the linear part of a Hamiltonian system of the form 

n 

^flincar = ^jPj 

i=i 

has no resonances of the form (|2.8p and the full system has no additional conserved quantities, as 
discussed in section 12.41 But in the near-resonance gives rise, at small nonlinearities, to the nearly- 
conserved quantity h, which allows for the dimension- reduction via averaging. 

Since the resonance above is not exact, the additional conservation laws are only approximate, and 
the the quantity h is not precisely conserved. Formally, one may perform a countable sequence of 
changes of variables that transform the system into a form that is completely integrable. In the limit, 
this corresponds to defining a change of variables given as a power series in e. Generally, this power 
series has radius of convergence zero, because the full system is not itself integrable, which we can see 
from the chaotic dynamics in the numerical solution given in figure [6?3h . This analysis suggests that 
the solution 16.3b is also very weakly chaotic, but with a much smaller chaotic region and a longer 
chaotic timescale. 

We are interested in the stability of the trivial solution of system (|7.ip : (pi^ps) = (0,0). This 
initial condition lies on the level set ft, = in iJroduccd- Thus, solutions to equation (|7.ip whose initial 
conditions satisfy ft 7^ cannot approach the origin and it suffices to set ft = in system (17. Sp when 
studying the stability of the trivial solution. Any stable or unstable manifolds to the origin must 
also lie in this level set. We observe what appears to be a near-homoclinic orbit in the numerical 
experiments presented in figure [^31 most clearly in row B, and by the above reasoning, any homoclinic 
orbit to the reduced system must be on the set ft = 0. Looking at the level set Hi{J, ; 0) =0 gives 
the following algebraic equation for the level set containing the origin. 

71^ + 72^^ +73(2 - J) cos-!/; = 0. 

This has the trivial solution J = as well as those that satisfy 

/ 71 +72 J .^ 
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The origin will have stable and unstable manifolds if this equation has a solution with J = 0. After 
some algebra, this simplifies to 



cosV' = 



which may only happen if 



Thus, there exist bifurcations at 



71 

73 



(ail22 — ^2222 + 02233) 1^ — S 



< 1. 



71 




(ail22 — ^2222 + a2233) — S 


73 




01223 



±ai 



223 



ail22 ~ ^2222 + 02233 



(7.5) 



which is simply a recapitulation of the bifurcation condition found by another method in equation (I5.3p . 
More simply, there exist fixed points of system (17.31) with sin?/; — and 



J tight 



2 (01122 ± 01223 ~ 02222 + 02233 ^ s/ 1^) 



— oiiii — 8aii22 — 4aii33 ± 4ai223 + 4a2222 — 802233 + 03333 



(7.6) 



Since J > by definition, these fixed points bifurcate from the origin exactly when the numerator 
vanishes, i.e. when the coefficients satisfy condition (|7.5p . For unstable values of v, the origin is always 
a multiple root, and thus is a nonhyperbolic fixed point of the averaged system. When ft, ^ 0, the 
origin is no longer a fixed point, but the fixed point defined by equation (|7.6p persists, for v > iyiiii{h), 
and this critical value now depends on h. For v below this bifurcation, the dynamics is described 
by monotonically decreasing angle 1/' (determined from the sign of the numerically calculated 71 ) and 
oscillating amplitude J. For v > i^hh, there exists a new fixed point, surrounded by a family of periodic 
orbits for which ip oscillates. This region is separated from the region of monotonic V' by a heteroclinic 
orbit connecting the line J = to itself. This difference can be seen by comparing parts (a) and (b) 
of figure 17.11 



Additional structure 

In addition, we note that the set Aovcn = {{J = ^,4')} is invariant under system (|7.3p when h = 0. 
Additional fixed points exist where 

cos^p^2l^tJl_ (7.7) 

73 

This corresponds to two fixed points a± on Aovcn- As in the case of the heteroclinic orbit given 
by (|7.4p . these exist only if the right hand side has magnitude less than one. This, then, gives a 
necessary condition on the amplitude v > i'f, for their existence, where 

lyp = ^ . (7.8) 

4O1111 — 01122 + 01133 =F 01223 — 02233 + 4O3333 

For v > vp, equation (|7.7p will have two solutions of saddle type connected by three heteroclinic 
orbits-two of them contained in Aevcn^and an additional fixed point Jicft of elliptic type with ijj — n 
and J = Jicft near i. The result of this bifurcation can be seen by comparing parts (b) and (c) of 
figure 17.11 Note that the left boundary J — corresponds to solutions on the odd invariant subspace 
ci = C3 = of system (|4.6I) . while the right boundary Aovcn represents solutions in the even subspace 
C2 = 0, and this figure shows a clear symmetry: above A'hh, there exists in the averaged equation an 
orbit homoclinic orbit to the odd subspace, and above N = evp, the averaged equations possess a pair 
of periodic orbits on the even subspace which are connected by three heteroclinic orbits. Note that 
the orbits a± correspond to periodic orbits of equation (|4.14p and cannot be found by the methods of 
section 14.21 
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Figure 7.1: The averaged {J,ip) phase plane with h = 0, and varying vahies of N. In (a), the phase 
changes monotonically. In (b), a new fixed point arises due to the HH bifurcation, and with it an orbit 
homochnic to the origin (J — 0). In (c), three new fixed points bifurcate from the fine J — ^. The 
colors represent the level sets of -ff average- 



When < h < 1, the two fixed points on the boundary of the phase disk persist for h small, with 
the bifurcation condition (|7.8p generalizing to 



i^Fih) = 1= (7.9) 

j{l + h)aiiii - (1 + h)aii22 + 01133 T V 1 - h?ai223. - (1 - ^)a2233 + 4(1 - ^)a3333 

Note that the averaged system (|7.3p is valid for small values of e and describes the dynamics for 
small |iV| demonstrated numerically in figure W% More concretely, the averaged system possesses the 
bifurcations described by equation (|5.3I) , but not the other two roots of equation (|5.2I) that may exist 
for N = (1). 

We end with a numerical computation that compares the integrable averaged system with the full 
system that displays Hamiltonian chaos. For the parameter values corresponding to figure [673b . column 
2, the averaged system has a phase space structure as in figure FTTb . Returning to Uj coordinates, the 
averaged system has four fixed points: the (Tright — vXight and CTicft — — \J Jioft are elliptic, and 
the two points a± on the boundary are hyperbolic. In the full system, h is not conserved, but the 
solution is still confined to a disk. In figure [7?2l we show, in blue, the Poincare section of figure [673h . 
along with several other solutions with the same value of N and H. On the left, there is a family of 
regular orbits surrounding cricft — with these parameter values, there is no chaos near this fixed points. 
Around the fixed point CTj-jght , we see what appears to be typical KAM breakup into Poincare-Birkhoff 
islands interspersed with the preserved KAM tori. As the parameter N increases, more and more 
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of the quasiperiodic orbits are destroyed, producing the fully-developed chaose, seen in figure 16.3d . 
The chaotic dynamics near the HH bifurcation have been analyzed recently in [51] • The dynamics 
produced by the HH bifurcation violate a "twist" condition assumed by the KAM theorem, so the 
structure of the system is not exactly the same as in the usual KAM setup. 




Reag ■ • • ■ 

Figure 7.2: Poincare map showing the chaotic region of figure 16.3b (blue), as well as many other 
features of the dynamics. The fixed points fJioft, right are given by yellow stars and a period-9 orbit 
inside the chaotic region in pink stars. The presence of a separatrix connecting cr± is clear on the left, 
and, on careful inspection, there is evidence for several other families of periodic orbits of period 3, 4, 
7, and 13. 



Comparison with other approaches 

Lahiri and Roy [43] considered the case of the non-semisimple HH bifurcation, so that the results they 
cite are not directly applicable. Using formal averaging of a different type, they find two types of 
bifurcations, depending on certain coefficients in the cubic and quartic terms in the Hamiltonian. In 
their Type 1 bifurcation, they find that there exists, for v > j'hh, a new fixed point of the averaged 
equations a distance d oc \/v — z^jh from the origin. Converting the solution (|7.6p to Cartesian 
coordinates, this is exactly what we find. 

In their Type II bifurcation, there exists no nonzero periodic orbit near zero on the unstable side of 
the bifurcation. This is the case for the bifurcations with A*" = O (1) that takes place between figure [673t 
rows (d) and (e) . Johansson makes the same observation in his study of the NLS trimer [10] but does 
not comment on the difference between the semisimple and non-semisimple bifurcations. 

8 Discussion and Conclusions 

While there have been a large number of papers, discussed in the introduction, examining the HH 
bifurcation and the onset of oscillatory instabilities in nonlinear wave equations, we believe this is the 
first to try to analyze the nonlinear dynamics that arise in such a system. In so doing, we discovered 
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the homoclinic orbit structure in the averaged equations, as weh as a new family of relative periodic 
solutions to the full system that exists when v > v-pili) in equation (j7.9p . 

While the analysis in the present paper is purely formal, we believe that the pieces are in place to 
make rigorous this paper's conclusions, as has been done in fT. In that paper, for a two- mode potential, 
it is shown that under suitable assumptions, that below some critical amplitude A/sb — Olf); there 
exists a steady solution with even symmetry which is stable to perturbations, and that for M > A/sb, 
there exist two new branches of asymmetric solutions which are orbitally Lyapunov stable. In ^6^, it 
is further shown that the time-dependent dynamics of the PDE solution are well-modeled, for long 
but finite times, by a finite-dimensional system that is equivalent to that of a particle moving in a 
potential which has just one well when J\f < A/sb but two wells when Af > A/sb- 

In the finite dimensional system of approximate equations near a symmetry-breaking bifurcation, 
derived in [H |6] , it takes one line of algebra to show the existence of the two new asymmetric solutions 
that are born when the bifurcation occurs. In system (|4.6p . the analysis is not so simple, and the 
new solution arising from the bifurcation appears, in its simplest form, as the fixed point (j7.6l) of 
Hamiltonian system ()7.3p that corresponds to a periodic orbit of system with Hamiltonian ()4.1ip . 
Proving the existence of this periodic orbit is a straightforward application of a paper from the late 
1980's by Chow and Kim 39 and will constitute the first step of a planned program to put the results of 
the present paper on a more rigorous footing. The chaotic dynamics near the HH bifurcation in a finite- 
dimensional system are rigorously demonstrated in [41!, 1421 . An attempt to rigorously demonstrate 
complex dynamics in NLS (jl.ip must start with an understanding how these results apply to the finite 
dimensional model (|4.14p . 

It should be noted that while in this system, it is possible to observe Hamiltonian chaotic motion, 
the underlying dynamics, given by system (|4.11l) are essentially two degree-of-freedom. Motion of 
such a system occurs on level sets of the Hamiltonian H which are three-dimensional manifolds in 
the four-dimensional phase space. Invariant tori in this system are two-dimensional subsets of these 
manifolds. The KAM theorem (or something very similar, see ^l]) implies that most of these tori 
persist when A/" — A/hh is small and positive. A two-dimensional torus separates the three-dimensional 
manifold, so that trajectories cannot cross from one side of the torus to the other. This implies that 
solutions starting near the odd-symmetric relative fixed points must remain near that point. If the 
linear system (j2.1l) is assumed to support a fourth eigenmode, with similar assumptions on the spacing 
of the eigenvalues, then in this weakly unstable regime, with six-dimensional phase space, solutions no 
longer need stay close to the fixed point, a process known as Arnol'd diffusion [36]. Further studies 
are planned to investigate this possibility. 

We have assumed throughout this paper that the potential V{x) enjoys even spatial symmetry. 



The HH bifurcation phenomenon discussed in this paper depends only on assumptions (|(Al)[)-(|(A4 



and not on this symmetry. Lacking such a symmetry, the finite-dimensional model (|4.1ip and its 
relative equilibria given in section [4.21 would be significantly more complicated, and the normal form 
for the HH bifurcation might no longer be semisimple. An interesting question would be to see how 
the dynamics change in the face of such asymmetry. 

Finally, when considered as a model for an optical waveguide, the system studied here should be 
straightforward to implement in a laboratory setting. Discussions are underway to make this happen 
and will form the basis of an experimental line of research. 
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A Some 3-soliton formulas 

For any three real numbers Kj satisfying ki > ^2 > K3 > 0, there exists a threc-sohton potential given 
by: 

u{x)^Mu{x)/V{x) 

and the three modes are given by 

V',(x)=AG(x)/I?(x) 

where 

— 2(k^ — K2)(ki — '«3)('^2 — K3)^(k2 + ^3)^ cosh2Kia; 
— 2(k^ — K2){k2 — Kf){Ki — K3)^(ki + K^)'^ cosh2K2a: 
— 2(k^ — K3)(k2 — K3)('«i ^ K2)^(ki + K2)^ cosh 2K3a:; 
+ K2)^(ki - K3)(k2 - K3)k3(ki + K3)(k2 + K3) cosh2(Ki - ^2)2; 
— (k1 — K2)'«2(ki + K2)(k2 — K3)(ki + Kz)^ {k2 + K3) C0Sh2(Ki — ^3)2; 

- K2)(k1 + '«2)(k1 - '«3)('«1 + K3)('«2 + ^3)^ COsh2(K2 " ^3)0; 
-(k1 - K2)^(ki - K3)(k2 - K3)k3(ki + K3)(k2 + K3) C0sh2(Ki + K2)x 
-{ki - K2)k2(ki + '«2)(ki - K3)^(k2 " ^z){k.2 + '«3) C0sh2(Ki + Kz)x 
-K,l{ni - K2){H1 + K2)('«l - K3){H2 - H^f {^1 + K3) COSh2(K2 + ^3)2;, 

Hi (.t) = (fi;2 + K3) cosh (k2 ~ '^3)2; + (^2 — K3) cosh (k2 + ^3)0;, 
H'2{x) = (ki + K3) sinh (ki — ^3)0; + (ki — K3) sinh (ki + K3)a;, 
N:i{x) — {ki — K2) cosh (ki + K2)x — {ki + K2) cosh (ki — K2)x, 

and 

X'(x) = (ki + K2)('«l + K3)(k2 - K3) cosh (ki - K2 - ^3)2: 
+ (ki - K2)(ki + '«3)(k2 + K3) cosh (ki + K2 - ^^3)^; 
+ (ki + K2)(ki - K3)(k2 + ^^3) cosh (ki - K2 + K3)x 
+ (k1 - '«2)(ki - K3)(k2 - '«3) cosh (ki + K2 + K3)x 

The three discrete eigenvalues are given by ^Ij = —k'j. 

B Remainder terms 

Equation (|4.4p depends on remainder terms i?2, -R3, and i?cont which we define here. The remainder 
terms in equations (|4.4ap - (|4.4cp are given by 

where Hj is given in equation (I4.2p and 

F = |ci*i +C2*2 +C3*3 +?7l^ (ci*l+C2*2+C3*3+f7)-|ci«'l + C2 5'2 + Cs^'s]^ (ci *i + C2 *2 + €3*3) 

The remainder term for the ri{x,t) equation (j4.4dp is given by 

-Rcont = —A/" • IIcontG 

where Ilcont is given in (|4.3p and 

G = |ci*i + C2*2 + C3*3 + + C2*2 + C3*3 + 
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